*******************************************************************************
****						Appendix Table G.18	- REGRESSION OF TOTAL TOLL ****
****			ON TIME DIFFERENTIALS: MODELS WITHOUT CONSTANT			*******
*******************************************************************************
use ".\data\clean\I10W_laneuse_dataset_15nov14_wcensus", clear

merge m:1 date hour using ".\data\clean\HV_ML_reliab.dta", keep(1 3) nogen 
gen reliabilityML=dist/p20_speedML-dist/p50_speedML
gen reliabilityHV=dist/p20_speedHV-dist/p50_speedHV
gen reliability_diff=reliabilityML-reliabilityHV
replace reliability_diff=0 if reliability_diff<0.01
keep if reliability_diff~=.
la var reliability_diff "Reliability"

*WTP calculation
gen TT_dif_hr=dist/MLspeed-dist/ELspeed
gen WTP2=charged_toll/TT_dif_hr
drop if ELspeed==.
drop if TT_dif_hr==.
drop if holiday==1
la var TT_dif_hr "Time in Hours"
drop if dow==0|dow==6
keep if acct_type=="PRIVATE"&occupancy~="HOV-3"
keep if hour>4 & hour<9
forval p=1/5{
	g TT_dif_hr`p'=TT_dif_hr^`p'
	}

estimates clear

*********************************************************
****	Col 1	Linear Model
*********************************************************
reg charged_toll TT_dif_hr  reliability_diff if TT_dif_hr>0 , cluster(rt_id) nocons
est sto a1
qui sum charged_toll if e(sample)==1
local toll1=round(`r(mean)',.01)
estadd scalar toll=`toll1'
qui sum TT_dif_hr if e(sample)==1
lincom	TT_dif_hr*`r(mean)'
local wtp1=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp1'

*********************************************************
****	Col 2	Quadratic Model
*********************************************************
reg charged_toll c.TT_dif_hr##c.TT_dif_hr  reliability_diff if TT_dif_hr>0 , cluster(rt_id) nocons
margins , at(TT_dif_hr=(0(.05).5))
marginsplot, name(WTPfromTT2_nocons, replace) yline(0) title(WTP from Travel Time - Quadratic Model without Constant) ///
	xtitle(Time Savings (in hours)) recast(line) recastci(rarea) plotopts( lc(black) ///
	 lp(dash) ) ciopts(fc(gs8%40) lw(none)) ///
	ytitle("Predicted WTP from Travel Time Savings" "(in dollars)") ///
		  xlabel(0(.1).5,format(%03.1f)) xtick(0(.05).5)
graph export .\results\appendix\figs\fig14_quadratic.png, replace		  		  
reg charged_toll TT_dif_hr TT_dif_hr2 reliability_diff  if TT_dif_hr>0 , cluster(rt_id) nocons
est sto a2
qui sum reliability_diff if e(sample)==1
scalar A0=`r(mean)'*_b[reliability_diff]
scalar A1=_b[TT_dif_hr]
scalar A2=_b[TT_dif_hr2]
forval p=1/2{
	qui sum TT_dif_hr`p' if e(sample)==1
	local mean`p'=`r(mean)'
	}
qui sum charged_toll if e(sample)==1
local toll2=round(`r(mean)',.01)
estadd scalar toll=`toll2'	
lincom	TT_dif_hr*`mean1'+TT_dif_hr2*`mean2'
local wtp2=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp2'
mata									
	a0 = st_numscalar("A0")
	a1 = st_numscalar("A1")
	a2 = st_numscalar("A2")
	s = polyroots((a0,a1,a2))
	sReal = select(s, Im(s):==0)
	sReal = Re(sReal)
	sGeq0 = select(sReal, sReal[.,.]:>0)
	sLeq1 = select(sGeq0, sGeq0[.,.]:<1)
	sol = sLeq1
	st_numscalar("solution", sol)
	st_local("solution", strofreal(sol))
end
local root2=round(60*`solution',.1)
estadd scalar root=`root2'

*********************************************************
****	Col 3	Cubic Model
*********************************************************
reg charged_toll TT_dif_hr TT_dif_hr2 TT_dif_hr3  reliability_diff  if TT_dif_hr>0 , cluster(rt_id) nocons
est sto a3
qui sum charged_toll if e(sample)==1
local toll3=round(`r(mean)',.01)
estadd scalar toll=`toll3'
forval p=1/3{
	qui sum TT_dif_hr`p' if e(sample)==1
	local mean`p'=`r(mean)'
	}
lincom	TT_dif_hr*`mean1'+TT_dif_hr2*`mean2'+ TT_dif_hr3*`mean3'
local wtp3=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp3'

*********************************************************
****	Col 4	Quartic Model
*********************************************************
reg charged_toll c.TT_dif_hr##c.TT_dif_hr##c.TT_dif_hr##c.TT_dif_hr  reliability_diff if TT_dif_hr>0 , cluster(rt_id) nocons
margins , at(TT_dif_hr=(0(.05).5))
marginsplot, name(WTPfromTT4_nocons, replace) yline(0) ///
	title(WTP from Travel Time - Quartic Model without Constant) ///
	xtitle(Time Savings (in hours)) recast(line) recastci(rarea) plotopts( lc(black) ///
	 lp(dash) ) ciopts(fc(gs8%40) lw(none)) ///
	ytitle("Predicted WTP from Travel Time Savings" "(in dollars)") ///
		  xlabel(0(.1).5,format(%03.1f)) xtick(0(.05).5)
graph export .\results\appendix\figs\fig14_quartic.png, replace		  		  		  
reg charged_toll TT_dif_hr TT_dif_hr2 TT_dif_hr3 TT_dif_hr4 reliability_diff  if TT_dif_hr>0 , cluster(rt_id) nocons
est sto a4
qui sum reliability_diff if e(sample)==1
scalar A0=`r(mean)'*_b[reliability_diff]
scalar A1=_b[TT_dif_hr]
scalar A2=_b[TT_dif_hr2]
scalar A3=_b[TT_dif_hr3]
scalar A4=_b[TT_dif_hr4]
forval p=1/4{
	qui sum TT_dif_hr`p' if e(sample)==1
	local mean`p'=`r(mean)'
	}
qui sum charged_toll if e(sample)==1
local toll4=round(`r(mean)',.01)
estadd scalar toll=`toll4'	
lincom	TT_dif_hr*`mean1'+TT_dif_hr2*`mean2'+TT_dif_hr3*`mean3'+TT_dif_hr4*`mean4'
local wtp4=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp4'
mata									
	a0 = st_numscalar("A0")
	a1 = st_numscalar("A1")
	a2 = st_numscalar("A2")
	a3 = st_numscalar("A3")
	a4 = st_numscalar("A4")
	a5 = st_numscalar("A5")
	s = polyroots((a0,a1,a2,a3,a4,a5))
	sReal = select(s, Im(s):==0)
	sReal = Re(sReal)
	sGeq0 = select(sReal, sReal[.,.]:>0)
	sLeq1 = select(sGeq0, sGeq0[.,.]:<1)
	sol = sLeq1
	st_numscalar("solution", sol)
	st_local("solution", strofreal(sol))
end
local root4=round(60*`solution',.1)
estadd scalar root=`root4'

*********************************************************
****	Col 5	Quintic Model
*********************************************************

reg charged_toll TT_dif_hr TT_dif_hr2 TT_dif_hr3 TT_dif_hr4 TT_dif_hr5  reliability_diff if TT_dif_hr>0 , cluster(rt_id) nocons
est sto a5
qui sum charged_toll if e(sample)==1
local toll5=round(`r(mean)',.01)
estadd scalar toll=`toll5'
forval p=1/5{
	qui sum TT_dif_hr`p' if e(sample)==1
	local mean`p'=`r(mean)'
	}
lincom	TT_dif_hr*`mean1'+TT_dif_hr2*`mean2'+ TT_dif_hr3*`mean3'+ TT_dif_hr4*`mean4'+ TT_dif_hr5*`mean5'
local wtp5=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp5'



esttab a1 a2 a3 a4 a5  using ".\results\appendix\tabs\ATC14.csv", replace  ///
nonumbers nodepvars ///
	cells(b(star fmt(%9.2f)) se(par)) star(* 0.10 ** 0.05 *** 0.01) ///
		stats(N r2  aic bic ll wtp toll root, ///
		fmt(%9.0f %9.2f %9.0f %9.0f %9.0f %9.2f %9.2f %9.4f) ///
		labels(Observations R-Squared AIC BIC Log-Likelihood "WTP from Travel Time" "Total Tol" ///
		 "Time Savings in Min. for WTP=0")) ///
	label mtitle("I" "II" "III" "IV" "V")
insheet using ".\results\appendix\tabs\ATG18.csv", comma clear
export excel using ".\results\appendix\AppendixFigsTabs.xlsx", ///
	sheet("Appendix Table G.18", replace)
putexcel set ".\results\appendix\AppendixFigsTabs.xlsx", sheet("Appendix Table G.18") modify
putexcel W1 = picture(.\results\appendix\figs\fig14_quadratic.png)
putexcel W50 = picture(.\results\appendix\figs\fig14_quartic.png)
